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Abstract. In this paper we propose a malaria within- host model with k classes of age for the 
parasitized red blood cells and n strains for the parasite. We provide a global analysis for this model. 
A competitive exclusion principle holds. If TZq, the basic reproduction number, satisfies TZg < 1, 
then the disease-free equilibrium is globally asymptotically stable. On the contrary if TZo > 1, then 
' generically there is a unique endemic equilibrium which corresponds to the endemic stabilization of 

, the most virulent parasite strain and to the extinction of all the other parasites strains. We prove 

, that this equilibrium is globally asymptotically stable on the positive orthant if a mild sufficient 

$Zh ' condition is satisfied. 
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1. Introduction. In this paper we consider intrahost models for malaria. These 
models describe the interaction of a parasite, namely a protozoa Plasmodium falci- 
parum, with its target cells, the red blood cells (RBC). During the past decade there 
has been considerable work on the mathematical modeling of Plasmodium falciparum 
infection [2, 14, 21, 22, 24, 23, 25, 28, 30, 52, 55, 56, 58, 64]. A review has been done 
J> ' by Molineaux and Dietz in [59] . 

(N : We give a brief review of the features of malaria. Malaria in a human begins 

' with an inoculum of Plasmodium parasites (sporozoites) from a female Anopheles 

mosquito. The sporozoites enter the liver within minutes. After a period of asexual 
reproduction in the liver the parasites (merozoites) are released in the bloodstream 
where the asexual erythrocyte cycle begins. The merozoites enter RBC, grow, and re- 
produce over a period of approximately 48 hours after which the erythrocyte ruptures 
releasing 8-32 "merozoites" daughter parasites that quickly invade a fresh erythro- 
cyte to renew the cycle. This blood cycle can be repeated many times, in the course 
of which some of the merozoites instead develop in the sexual form of the parasites: 
gametocytes. Gametocytes are benign for the host and are waiting for the mosquitoes. 

The first mathematical model of the erythrocyte cycle was proposed by Anderson, 
May, and Gupta [3]. This original model has been extended in different directions 
[2, 3, 21, 25, 28, 30, 64]. 

The original model [3] is given by the following system: 

X = A — /iajX — (3x m, 
(1.1) <( y ^ I3xm- fiyy, 

m — r fly y — fim m — f3 x m. 
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The state variables are denoted by x, y, and m. The variable x denotes the concen- 
tration of uninfected RBC, y the concentration of parasitized red blood cells (PRBC), 
and m the concentration of the free merozoites in the blood. 

We briefly sketch the interpretation of the parameters. Parameters iJ,x, l^y, and 
lim are the death rates of the RBC, PRBC, and free merozoites, respectively. The 
parameter /3 is the contact rate between RBC and merozoites. Uninfected blood cells 
are recruited at a constant rate A from the bone marrow and have a natural life- 
expectancy of — days. Death of a PRBC results in the release of an average number 
of r merozoites. Free merozoites die or successfully invade a RBC. 

This system is isomorphic to numerous systems considered in the mathematical 
modeling of virus dynamics; see [60, 61, 62] and the references therein. Some authors 
ignore the loss term —pxm that should appear in the m equation. Indeed without 
this loss term, merozoites can infect RBC without themselves being absorbed, and 
this allows one merozoite to infect more than one RBC. 

The original and the derived malaria models were intended to explain observa- 
tions, namely parasitaemia, i.e., the concentration y of PRBC and also the decrease of 
the healthy RBC leading to anaemia. An important characteristic of Plasmodium fal- 
ciparum, the most virulent malaria parasite, is sequestration. At the halfway point of 
parasite development, the infected erythrocyte leaves the circulating peripheral blood 
and binds to the endothelium in the mierovasculatiire of various organs where the cy- 
cle is completed. A measurement of Plasmodium falciparum parasitaemia taken from 
a blood smear therefore samples young parasites only. Physician treating malaria use 
the number of parasites in peripheral blood smears as a measure of infection, and 
this does not give the total parasite burden of the patient. In some respects this is 
a weak point of the model (1.1). Moreover antimalarial drugs are known to act pref- 
erentially on different stages of parasite development. These facts lead some authors 
to give a general approach to modeling the age structure of Plasmodium parasites 
[22, 23, 24, 57]. Their model is a linear catenary compartmental model. This model is 
based on a finite number of compartments, each representing a stage of development 
of the parasite inside the PRBC. The models describe only the dynamics of the mor- 
phological stage evolution of the parasites and make no allowance for the dynamics 
of the healthy RBC. 

In this paper we propose a model which combines the advantages of the two 
approaches. We also consider this model with different strains for the parasites. To 

encompass the different models of the literature we allow, in this model, to ignore or 
not the loss term in the m equation. To begin we consider the model with one strain: 

X = f{x) - - jSxm, 
yi = /3xm- aiyi, 

m = 7iyi - "2 2/2, 

ilk =7fe-i2/fe-i -oikyk, 
rn = r^kyk- l^mm-uPxm. 

In this system f{x) — n^x is the density-dependent growth rate of RBC. The other 
parameters arc positive. In the model of Gravenor et al. [21] ai = 7i -f /Ltj, and hence 
ai > 7i. We do not need this requirement, which implies that our model is not 
necessarily a catenary compartmental model. In the literature the parameter u takes 
the values u = when the loss of the merozoite when it enters a RBC is ignored or 
takes u = 1 when this loss is not ignored. In our analysis u is simply a nonnegative 



(1.2) 
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parameter. Except for these generalizations this system has already been suggested 
by Gravcnor and Lloyd [21] in their reply to the criticism of Saul [64]. We provide 
a global analysis of this system related to the basic reproduction ratio TZq of the 
considered model. 

One problem is how to decide upon the number of parasite compartments in 
the model. A starting point can be the morphological appearance of the parasite. 
But if the objective is to reflect the distribution of cycle lengths, the number of 
compartment can be increased to obtain a gamma distribution. Finally the two 
approaches can be combined: some compartments are for morphological reasons and 
others are for behavioral reasons. Then this model can also be interpreted as the 
application of the method of stages (or the linear chain trick) to the life cycle of 
PRBC [3, 31, 47, 49, 48, 51]. In other words a chain of compartments is included to 
generate a distribution of lags. It is also possible to add a class tjk+i in order to allow 
for the production of gametocytes. Different numbers of stages, ranging from 5 to 48, 
are used in [20, 22, 23, 24]. 

It is well grounded that a falciparum infection consists of distinct parasite geno- 
types. The model of Anderson, May, and Gupta has been extended in this direction 
[25, 66] . With regard to such features we propose a model with k stages for the infected 
RBC, production of gametocytes, and n genotypes, in the population of parasites. 

One of the important principles of theoretical ecology is the competitive exclusion 
principle which states that no two species can indefinitely occupy the same ecological 
niche [7, 8, 11, 17, 25, 39, 53, 54]. We provide a global analysis of this model and 
obtain a generic competitive exclusion result within one host individual. This confirms 
the simulation results obtained in [25]. We compute the basic reproduction ratio TZq 
of the model. For this model there is always a disease- free equilibrium (DFE). To put 
it more precisely this equilibrium corresponds to the extinction of all the parasites, 
including the free parasites and the intraerythrocyte parasites. We prove that if 
7?.o < 1, then the DFE is globally asymptotically stable (GAS); in other words the 
parasites are cleared. If TZq > 1, then, generically, a unique endemic equilibrium exists 
corresponding to the extinction of all the strains of parasites but one. We prove that 
this equilibrium is GAS on the positive orthant under a mild condition. For example 
this condition is automatically satisfied when u = and f{x) = A — /ij; .t. When u ^ 
the criteria, obtained for deciding the winning strain, differs from other results in the 
literature. To each i-strain can be associated a basic reproduction number TZq and a 
threshold To*- It turns out, when u 7^ 0, that this is precisely this threshold Tq which 
distinguishes the fate of the strain and not TZ^ at the difference of [7, 11]. 

The paper is organized as follows. In section 2 we introduce the model with k 
stages for the infected RBC and one parasite strain, with and without gametocyte pro- 
duction. We compute the basic reproduction number and provide a stability analysis. 

In section 3 we consider the model of Anderson, May, and Gupta with n distinct 
genotypes and production of gametocytes. This model with a constant recruitment 
function for the erythrocytes, two strains, and one class of age has been proposed 
in [25]. We have studied this model in [1]. Here using the computation of section 
2, we prove for the general n strain k class of age model that if TZq < 1, then the 
parasites are cleared and ii TZq > 1, then generically the different genotypes cannot 
coexist. Namely a unique equilibrium exists, for which only one genotype is positive, 
and which is GAS on a dense subset of the nonnegative orthant. This result confirms 
the simulations given in [25]. 

Global results of stability for the DFE as well for the endemic equilibrium for 
epidemic models are not so common [26, 27, 33, 43, 65, 67, 68]. Global stability 
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results for the endemic equilibrium have often been obtained by using monotone 
system techniques [29, 36]. Usually the Poincare-Bendixson property of monotone 
systems in dimension 3 is used [40, 41, 42, 43, 44, 45]. Our results generalize the 
results of [13]. 

2. Stability analysis of a one strain model with k stages. We consider 

a general class of systems. The haemopoiesis is a complex system. In the cited 
references the recruitment of RBC is given hy A — x. In this paper we will use a 
more general function ^p{x). In a more complex system the haemopoiesis could be an 
input coming from another system: 

X = f{x) — fjbxX — j3xm = ip{x) — I3xm, 

yi = /Sxm- ai yi, 
m = 71 2/1 - "2 2/2, 

ilk = 7fc-i Vk-i - cxk yk, 
m = rjkyk-IJ'mm-ul3xm. 

We denote by y the column vector {yi, . . . ,yk)'^- The parameter u is nonnegative. 
The reason for this parameter is to encompass some malaria models in which the 
term —fixm can appear or not. In [2] Anderson has considered a system without 
the —(3xm in the m equation. In [60] all the basic models of virus dynamics are also 
without this term. One feature of Plasmodium falciparum, responsible for the deadly 
case of malaria, is that more; than one parasite can invade RBC. In this case u is the 
mean number of parasites invading RBC and thus disappearing from the circulating 
blood. 

Some authors [25, 56] have included in the model production of gametocytes. In 
the course of the production of merozoites from bursting erythrocytes, some invad- 
ing merozoites develop into the sexual, nonreplicating transmission stages known as 
gametocytes. The gametocytes are benign and transmissible to mosquitoes. We can 
also, following these authors, include a production of gametocytes in our model. If 
we denote by yk+i the "concentration of gametocytes," the model becomes 

X = f{x) — HxX — (3xm = (p{x) — Pxm, 
yi = Pxm-a-iy-i, 
m = 71 2/1 - "2 2/2, 

Vk =7fe-i2/fc-i -akVk, 

Vk+i = Plkyk- afe+i J/fc+1, 
TO = r 7/c y/c - m - M /3 x m. 

We start to analyze the system with minimal hypothesis on / but nevertheless plausi- 
ble from the biological point of view. The function / gives the production of erythro- 
cytes from the bone marrow. The function </5(x) = f{x) — models the population 
dynamic of RBC in the absence of parasites. The RBC have a finite lifetime, and 
then represents the average per capita death rate of RBC. The function / models 
in some way homeostasis. In this paper we suppose that / depends only on x. It 
could be assumed that the recruitment function depends on x and the total popula- 
tion of erythrocytes x + ^^ yi. In this paper we will analyze the simplified case which 
is the model considered in all the referenced literature. The rationale behind this 
simplification is that in a malaria primo-infection typically y is in the order of 10~^ 



(2.1) 
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to 10~* of the concentration of healthy erythrocytes x. This can be confirmed from 
the data of malaria therapy. In the last century neurosyphilitic patients were given 
malaria therapy, which was routine care at that time. Some of them were infected 
with Plasmodium falciparum. Data were collected at the National Institutes of Health 
laboratories in Columbia, SC and Milledgeville, GA during the period 1940 to 1963 
[12]. 

We assume that / is a C^. Since homeostasis is maintained we assume that the 
dynamic without parasites is asymptotically stable. In other words, for the system 

X = f{x) - i_ixx = (p{x) 

there exists a unique x* > such that 

(2.3) <^(a;*) = 0, and (p{x) > for < a; < x*, and ip{x) < for a: > x*. 

2.1. Notation. We will rewrite systems (2.1) and (2.2) in a condensed simpler 
form. 

Before we introduce some classical notation. 

We identify vectors of M" with n x 1 column vectors. ( | ) denotes the euclidean 
inner product. \\z\\2 = {z \ z) is the usual euclidean norm. 

The family {ei, . . . , €„} denotes the canonical basis of the vector space M". For 
example ei = (1,0,..., 0)-^. We denote by the last vector of the canonical basis, 
e„ = (0,...,0,1)^. 

li z € R", we denote by Zi the ith component of z. Equivalently Zi — { z \ Ci) . 

For a matrix A we denote by A{i,j) the entry at the row i, column j. For matrices 
A, B we write A< B ii A{i, j) < B{i,j) for all i and j, A < B if A< B and A^ B, 
and A<^ B if A{i,j) < B{i,j) for all i and j. 

A^ denotes the transpose of A. Then {zi \ Zi) = zf Z2- The notation A~^ will 
denote the transpose of the inverse of A. 

For this section we rewrite the systems (2.1) and (2.2) under a unique form: 



(2.4) 



X = f{x) - l3x{eoj\ z), 

z = /3x {ea, I z) ei + Aqz - upx {cuj \ z) e^- 



In the case of the system (2.1) we have for Aq 



(2.5) 



—ai 

7i -"2 
72 








-as 










-ak 



and an analogous formula for (2.2). 

We define the matrix A{x) — Ao — (3 x e^j e^. This a Metzler stable matrix. (A 
Metzler matrix is a matrix with nonnegative off-diagonal entries [5, 32, 50].) 

It is not difficult to check that the nonnegative orthant is positively invariant by 
(2.4) and that there exists a compact absorbing set K for this system. An absorbing 
set D is a. neighborhood such that a trajectory of the system starting from any initial 
condition enters and remains in D for a sufficiently large time T. 



ANALYSIS OF NEW MALARIA INTRAHOST MODELS 



265 



2.2. Global stability results. We can now give the main result of this section. 
Theorem 2.1. We consider the system (2.4) with the hypothesis (2.3) on (p 
satisfied. We define the basic reproduction ratio of the system (2.1) and (2.2) by 



(2.6) 



Ho = 



71 • • • 7fc 



Hm + u (3 X* ai ■ ■ ■ ak 



1. The system (2.1) is GAS on ]R^+^ (respectively, (2.2) on R^+^j at the DFE 
(a;*, 0, . . . , 0) if and only ifTZo < 1- 

2. If TZo > 1, then the DFE is unstable and there exists a unique endemic 
equilibrium (EE) in the positive orthant, {x, z) » 0, given by 



(2.7) 



tlr. 



13 



71 • • • 7fe 



Q!i • • • afe 

^ z = ip{x) (-Ao)"^ (ei - ue„). 



Denoting a* = - max^-gjo, iv'ix) ), if 

(2.8) M,/3^(,x-) < aVm, 

then the EE is GAS on the nonnegative orthant, except for initial conditions on the 
X-axis. 

Proof of Theorem 2.1. To begin we will consider the system (2.1) without game- 
tocytes, i.e., the system (2.4) with Aq as defined in (2.5). The stability analysis for 
(2.2) follows easily from the stability analysis of (2.1). 

In a first step we will compute TZq. We use our preceding notation and define 
A* = A{x*), i.e., the matrix computed at the equilibrium x* of (p, which is a stable 
Metzler matrix. We will use, repeatedly in what follows, the property that if M is a 
stable Metzler matrix, then —M^^ > [5]. The expression of TZo is obtained easily 
by using the next generation matrix of the system (2.1) [9, 15, 16]. We have for the 
basic reproduction number 



ei 



nQ = 0x* ( - {A 



If we remark that the matrix A* is the matrix Aq modified by a rank-one matrix, 
namely A* = Aq — u /3 x* Cu, e"^, we can use the Sherman-Morrison- Woodbury formula 



uj3x* 



l + upx*el (-^o)"' 



(-^o)"' e^el (-Ao)-^ 



or equivalently 



ujix* 

IJ.m+ (3x* 



eu, el {-Ao) ^ . 



This shows that — (A*) ^ is obtained from — ^ multiplying the last line of — j4q ^ 

- (-4o)~^ ei I ea 



by ^^-Tr/ax- Then we get 



■Ro = px* 



u X* 



and then in computing the last entry of the first column of Aq we obtain (2.6). 
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We remark that 7?.o > 1 is equivalent to the following threshold condition: 

To = ( - (^o)~^ ei I \ - u = j3x* {Ao)~^ (ei - u e„) | e^, ) > 1. 

We are now ready to analyze the stability of the DFE. 

It is well known that if 7?.o > 1, then the DFE is unstable [15], which implies that 
the condition T^o < 1 is necessary for stability. 

To prove the sufficiency, in a second step, we consider the following function 
defined on the nonnegative orthant: 

(2.10) VnFEiz) = /3x*{e^\ {-Ao^)z). 
Its time derivative along the trajectories of system (2.4) is 

Vdfe = Px{euj \ z) ^x* {e^ \ (-^o)~^ {ei - ucu,)) - ^ x*{eu; \ z) 
or equivalently, using the expression of To given in (2.9), 

(2.11) VDFE = P{eo,\z} (Tox-x*). 

Now we take as a candidate Liapunov function, defined on the nonnegative orthant 
minus the hyperplane face x = 0, 



V = {x-x* Inx) - x*{l - Inx*) + Vdfe{z). 

This function is positive definite (relatively to the DFE) on IR^"t^^Q = {{x,y,m) e 
M^"^^ : a; > 0}. Its time derivative is given by 

V=^^^^x)^^x-x*)P{e^\z) + p{e^\z) {%x-x*) 

X 

or assuming TZq <1 

l> = ^^^(x) + /Jx(e„ |z) (To-l) < 0. 

X 

By assumption (2.3) we have (x — x*)(p{x) < for all a; > 0. Therefore F < for all 
{x,z) € R^^^^Q, which proves the stability of the DFE. Its attractivity follows from 
LaSalle's invariance principle [6, 37, 38], since the largest invariant set contained in 
{{x, z) e M+|cc^>o : V = 0} is reduced to the DFE. On the other hand the vector field 
is strictly entrant on the face a; = 0. Hence the whole orthant M^"^^ belongs to the 
region of attraction of the DFE. 

Now we assume that '1Zq> 1. The equilibria {x,z) of the system, different from 
the DFE, are determined by the relations 

z = px{z\e^) {-Ao)~'^ {ei-ue^). 

Replacing z in ( z | ) we obtain 

(2.12) {z\e^)=px{z\e^){{-Ao)-\ei-ue^) \ e^). 
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If I e^^) = 0, then (p{x) = 0, we obtain x = x*, and hence ^ = 0; i.e., the 
corresponding equilibrium is the DFE. In the other case, i.e., {z \ e^^) ^0, the 
relation (2.12) gives 



(2.13) 

Using ((-Ao)"^ Cu 



X ( (-^o) ^ (ei -ueu)\euj) = 1. 



we finally have 



x" 
To' 



13 ^J,m ( (-^o)~^ ei I ) 
We deduce that if TZq > 1, then < a; < a;*, and hence (fi{x) > 0. Therefore 

z = ip{x) {-Ao)~^ (ei - ueu). 
The last component of ^, ( ^ | 6;^ ) = m, is given by 

(p{x) 



m 



px 



> 0. 



The k first components of z are given by the k first components oiip{x) {—Aq)~^ ei. It 
is straightforward to check that the first column of (— ^o)~^ namely (— ^o)~^ Ci ^ 0, 
which proves that z ^ 0. We have then proved that there is a unique EE in the 
positive orthant if and only if T^o > 1- 

Finally we will prove a sufficient condition for the global asymptotic stability 
of the EE. To this end we define the following candidate Liapunov function on the 
positive orthant minus the face corresponding to a; = 0: 

k 

(2.14) VEE{x,y,m) = a{x - xlnx) + ^ {yi - yjlnj/j) + bk+i {m-fh Inm). 



This function has a unique global minimum in {x, y, fh) . We will choose the coefficients 
a,bi,bk+i such that in the computation of V, the linear terms in yi and m and the 
bilinear terms in xm cancel. Let us show that it is possible with positive coefficients. 
To this end we rewrite the function Vee using the notation z = {y,m)'^, In^ = 
(In2;i,ln2;2, . . . ,ln^fc+i)^, and b= (6i, . . . , 6fe, ^fe+i)"^: 

Vee{x,z) = a{x — x\a.x) + {b\ z — diag{z)hi z) . 

Consider the block matrix 



M 



-1 

Pxe^ 



^0 



Using classical Schur complement techniques and the relation (2.13) on x, we have 

det(M) = det(Ao)[-l + ^5(ei - uca,)"^ (-^0 ^) e„] 

= det(Ao)[-l + /35(-Ao^(ei-Me„) | ej)] =0. 

Since the matrix Ad is obviously of codimension 1 (Ao is nonsingular) the kernel of 
M is of dimension 1. Then there exists a G M and b € such that 



(2.15a) 



a = (ei - u Ca,)^ 6 = ( 6 I ei - u e^, ) 
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and 
(2.15b) 



Since the kernel is one dimensional, a can be chosen arbitrarily. Thanks to the struc- 
ture of Aq, if a > 0, then 6 0. 

The derivative of V along the trajectories of (2.4) is given by 



Vee = a '■ 



■ tp{x) - a(3x{ei^ \ z) + a {3 x {e^j \ z) + j3x{eu \ z) {b \ ei - ue^j) 



+ {b\ Aqz) + {b\ diag(z) diag(2;)-^ z ) 

^ ^ 

= a (p(x) + (b\ diag(z) diag(z)~"^ z ) 

X 

+ aj3x {e^ \ z) + {b \ Aqz ) + f3x{e^^ \ z) [{b \ e-i - ue^) - a) . 

Using the relation (2.15b) we see that 

{b \ Aqz) = -a/3s((^o'^)eu, | Aqz) = -apx{e^ \ z). 

Therefore the linear terms in z cancel. The same is true for the bilinear terms thanks 
to the relation (2.15a). Finally we get 



EE = a ■ 



x 



■ ip{x) + {b\ diag(2:) diag(z) ^ i ). 



We choose bk+i = 1 = {b \ e^) = al3x{—A e^^ \ e^) = a[}x-i^. In other words 
a = With the hypothesis T^o > 1 we have a > 0, and hence 6 as wanted. 
With this choice developing V gives 



Vee =af(x) - afx^x- af{x) - + a /j.^ x - bi Pyi bai-iyi-i — 

X Vi ~^ Vi 



Em 
Oi OLi yi - r^k Vk Vujimx + 
m 



We collect some useful relations between our coefficients at the EE. We have from the 
definitions of a and b, since bk+i — 1, 



(2.16) 



a + u = bi, 
bi ai = 7i 62, 
62 0:2 = 72 bs, 

bk-i otk-i = 7fe-i bk, 
bkOik = r^k- 



Prom these relations and the properties of the EE z we have 

(2.17) bil3xfh = biaiyi = bi-ji-iyi-i =rjkyk 
and 

(2.18) aaiyi = Hmm. 
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Replacing, in the expression of V, a^^^ by a f{x) — aPxm = a.f{x) — aai yi we obtain 

X X 

Vee = krjkVk + a fix) + af{x) + [u^xifi - afi^x)- - af{x)- 

_ X m yi ^ _ y, _ yk m 

-bipxm-- > Oiji-iyi-i- r-fkyk- • 

X m yi ^ yi-i yi yk m 

1=2 

Using again the relations between the coefficients we get 

X X 

Vee = kr-fkVk + af{x) + af{x) + [r^kVk - af{x))- - af{x)- 

X X 



k 
i=2 



X m yi 
X fh yi 



Vi-i Vi - Vkfri 

nkVk- — 

Vi-i Vi Vk m 



and finally 



Vee = a 



f{x) + f{x)-f{x)^-f{x)^ 



■rjkVk 



Vi-i Vi Vk m 



X X m yi ^ yi-i yt yk m 



Now we will use the fact that there exists ^ in the open interval ^ e ]x,x[ such that 
f{x) = f{x) + {x — x) f'{0- Replacing in the preceding expression gives 



Vee = af{x) 



2- - - - 

X X 



-r^kVk 



X X m yi 



E 



Vi-i Vi Vk m 



X X m yi ^ yi-i yt yk m 



Using the relations (2.16)-(2.17) we have 

af{x) = (bi - u)f{x) = bi{na:X + fixfh) - uf{x) = bi fi^ x + rjkVk - uf{x). 
Replacing in the preceding expression of V gives 



Vee = {h iixX-uf{x)) 



2 - - - - 

X X. 



+ rjkyk 
This can also be written 



X X m yi 
X X ffi yi 



i=2 



Vi-i Vi _ m^rn 
Vi-i Vi Vk m 



Vee = ^{x, y,m) = - [bi HxX-u f{x) - ax f '{£,)] 



{x — a;)^ 



XX 



(2.19) 



-rjkUk 



X _ X m y]_ _ si-^ Vi-i yi _yk'm 
X X fh yi ^ yi-i yi yk m 
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The term between brackets in the last expression of V is nonpositive by the inequality 
between the arithmetical mean and the geometrical mean. Then a sufficient condition 
for y < is 

biii^x-ufix) -axf'i^) > 0. 

Moreover with this condition V is negative, except at the EE for the system (2.1). 
This proves the global asymptotic stability of the EE on the positive orthant for the 
system (2.1). 

The vector field associated with the system is strictly entrant on the faces of the 
orthant, except the x-axis, where it is tangent. The basin of attraction of the EE is 
then the orthant, except the x-axis, which is the stable manifold of the DFE. 

Using the function (p{x) = f{x) — the preceding condition is equivalent to 

uip(x) < —ax(p'{^), 
or equivalently, replacing a by its value ci = ^ > the condition becomes 

Setting a* = — maXa,g[o,x*] a sufficient condition for global asymptotic stability 

of the EE is 

TZo > 1 and u(}ip{x) < fimOt* ■ 

We have proved the theorem for the system without gametocytes. We have seen 
that TZq does not depend on the production of gametocytes. If 7?.o < 1, it is easy, 
integrating the linear stable yk+i equations of (2.2) from the solutions of (2.1), to 
see that the DFE is asymptotically stable and that all the trajectories converge to 
the equilibrium. The same argument is used when TZo > 1. This ends the proof of 
Theorem 2.1. □ 

Remark 1. If this model is a model for a within- host model of malaria, each 
coefficient ai is made of the mortality of the z-class and the rate of transmission 
in the i + 1-class: = /i^ + 7^. This implies that 7^ < ai. Wc do not need this 
assumption, and our conclusions are valid for our more general model. The only 
hypothesis is that the parameters of the system are positive. 

Remark 2. In the proof of Theorem 2.1 the quantity 

13 x* (Ao)'^ (ei - ue^) \ Cuj^ , 

which we have called To when TZo > 1, plays a prominent role. When TZo < 1 and 
three cases occur: G<7o<lor7o<0or7o = 0. 

In the two first cases we can define x = j^, and we obtain an equilibrium {x,z) 
of the system which is not in the nonnegative orthant (either i < or z < 0). 

In the third case, the computations, done in the proof of Theorem 2.1, for the 
research of an equilibrium show that {z \ e^) =0, and hence z = 0, and finally the 
equilibrium is the DFE (a;*,0). 

We introduce a definition of To that will simplify future computations. The case 
To = is special, since To = ^ is no longer true. However this case can be thought, 
by convention and misuse of language, as a; = +00. 
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Definition 2.2. We define for the system (2.1) the threshold 
(2.20) To = ^ =Px* (- {Ao)-' (ei -ue^)\ e^ 



^ 71 ■■■Ik 
tti • • • a/c 



When To ^ we have also To = ^ . 

Remark 3. It should be pointed out that the kind of Liapunov function defined 
by (2.14) has a long history of application to Lotka-Volterra models [18, 19] and was 
originally discovered by Volterra himself, although he did not use the vocabulary and 
the theory of Liapunov functions. Since epidemic models are "Lotka- Volterra" like 
models, the pertinence of this function is not surprising. Similar Liapunov functions 
have been used in epidemiology [4, 34, 35, 46, 63], although with different parameters. 
We have already used this kind of function in a simplified version of this paper in [1]. 

2.3. Comparison with known results. Our stability result improves the one 
of De Leenheer and Smith [13] in two directions: 

1. We introduce n stages for latent classes. 

2. Our sufficient condition for the global asymptotic stability of the endemic 
equilibrium is weaker than the one provided in [13]; for instance the sufficient condi- 
tion given in Theorem 2.1 is satisfied for malaria parameters given in [3], while the 
condition of [13] is not satisfied. 

2.4. Application to the original AMG model [3]. The original Anderson- 
May Guptka model is a three dimensional system (1.1) which has the same form as 
system (2.1) with f{x) = A. The sufiicient condition (2.8) applied to the AMG model 
(1.1) can be written 



(2.21) 



/3A < 



1 



fix l^n 



For the system (1.1), it is possible to give a weaker sufficient stability condition. 
Proposition 2.3. IfTZo>l and I3A < (-v/r + Vr- 1)^ /j-x jJ-m, then the EE is a 

GAS steady state for system (1.1) with respect to initial states not on the x-axis. 
Since in general the parameter r is larger than 2 (see, for instance, [28]), we have 

(v^ + v^)'>7^. 

Proof Thanks to the computations done before, we have for system (1.1) 

, ^ \^ X XI I X y 7h X m y 

VEE = ir-l)A 2---- +riiyy 1 + --^ - — - 

L X xi I X y m X m y 

Define X = ^ and 5 = ¥ — . Then one can write 

X y m 

Vee = -(r - 1)A^^^ +rtXyy{l + X-S-§) 

= -{r - 1)A^^ + r y *(X, S). 

We have *(A, S)>O^X<S<lorX>S>l. On the other hand *(X, S) < 
*(X,\/X) = {VX-lf. Therefore 



(2.22) 



Vee < (r - 1)A{VX - If + 1 + ^) - 1 - ;;^) • 
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We have ^ 



< 1. Hence for X < X* = ^ = ^^^^x* we have the 



following: 



1 



1 



< 



— 1 < 0, since by assumption 



— < {y^+^yr — 1)^. Therefore, the derivative of Veb along the trajectories 
of system (1.1) is negative definite on the set = {(a;, ?/,m) e : < a; < a;*, y> 
0, m > 0}. By continuity, there exists e > such that Vee is negative definite on 
the set T>e = {{x,y,m) € R\ : < x < x* + e, y > 0, m > 0}. The global 
asymptotic stability of the EE follows from the fact that is an absorbing set for 
system (1.1). □ 

3. The general case: n strains with k classes of parasitized erythro- 
cytes. We define the following system with k classes and n parasite strains: 



(3.1) 



X = fix) - fixX - X ^ A mi = (p{x) - x'^ j3i rrii 



i=l 

and for i = 1, . . . , n, 
yi,i = Pixrrii - an 

y2,i = 7i,j yi,i - Q!2,jj/2,j, 

yk,i = ik-i,i yk-i,i - oik,i yk,i, 

9i = kyk,i - fJ-gi 9i, 

rrii = Ti 7fe,i yk,i - jJ-rm rrii - uPixrrii. 



1=1 



As in preceding sections we rewrite the system as 



(3.2) 



X = (p{x) - x"^ I3i{zi\ ej,„ ) 
and for i = 1, . . . , n, 



where the matrix Ai is the analogous of the matrix Aq defined in section 2.2, but 
corresponding to the genotype i, and the vectors ej^i and Ci^^ are defined accordingly. 
We drop the index in >1 for readability. 

Theorem 3.1. We consider the system (3.1) with the hypotheses (2.3) satisfied. 
We define the basic reproduction ratio TZo of the system (3.1) by 



TiPiX* 



7i,» • • • 7fc,i 

Mm* +uPiX* ai,i---afc,i 



and 



TZo = max TZq. 



1. The system (3.1) is GAS on R_^_ at the DFE (a;*, 0, . . . , 0) if and only if 
■Ro<l. 

2. If TZo > 1; then the DFE is unstable. If Rq > 1, there exists an EE in the 
nonnegative orthant corresponding to the genotype i, the value for the other indexes 
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j ^ i are yj = rrij = 0, and 



(3.3) 



9i = 



* L • • • otk,i 
Si _ 



where we denote by Zi^k the fcth component of Zi . 

3. We assume TZq > 1. We define Tq as in Definition 2.2. We assume that the 
generic conditions Tq ^ Tq are satisfied for i ^ j- We suppose that the genotypes 
have been indexed such that 

To' >T^ >--->7^. 

Then the EE corresponding to xi is asymptotically stable and, the EEs corresponding 
to Xj for j 7^ 1 (for those which are in the nonnegative orthant) are unstable. 

4. We assume that the preceding hypothesis Tq > Tq is satisfied with TZq > 1. 
We denote it by a* = — maLXxe[o,x'] {f'i^) )• Then if 

the equilibrium (xi, t/i, mi, gi, 0, . . . , 0) is GAS on the orthant minus the x-axis and 
the faces of the orthant defined by yi = mi = gi = 0. In other words the most virulent 
strain is the winner and the other strains go extinct. 

Proof. As in Theorem 2.1 there exists a forward invariant compact absorbing set 
in the nonnegative orthant for the system (3.1), and hence all the forward trajectories 
are bounded. The variables gi do not affect the dynamical evolution of the variables 
X, yij, mi, and so we can consider the system without the production of gametocytes. 
We use the Liapunov function 



Vdfe{z) = ^ VoFEiZi) = |3^ X* { Bi,^ | {-A^ ')Zi ). 

i=l i=l 

Using the system written as (3.2) and the computation (2.11) we easily obtain 

n 

Vdfe = ^ Pi{ ei.o, \zi) (To a; - x*) . 

i=l 

Now we define the Liapunov function on the nonnegative orthant minus the hyper- 
plane face a; = 

n 

V{x, z) = {x — x* Inx) — x*{l — Inx*) + VoFEizi) 



which gives 



V = 



x — X 



^i^) + 51 ^*^' ^ I ^'''^ ) - xl3i{Zi\ Bi^oj ) 



i=l 



+ Yl Pi{^i^<^ I Zi){Vx-X*) 
1 

n 

'^{x) + ^ ft ( B^^u \Zi)x [Tq - 1) 



i=l 

a; — a;" 
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Since TZq < 1 for all index i, we have Tg < 1, and hence V < 0. The conclusion 
follows by Lasalle's invariance principle and consideration of the boundary of the 
positive orthant. 

Now we assume T^-o > 1. The instability of the DFE follows from the properties 
of TZq [15]. Wc assume that the genotypes are indexed such that their corresponding 
threshold are in decreasing order Tq > Tq > ■ ■ ■ > Tq- 

We will define a Liapunov function on the nonnegative orthant minus the manifold 
defined by the equations x = yi = rrii = 0. For this we need to recall the definition 
of the function VEE{x,yi,mi) defined in (2.14): 

k 

VEE{x,y,m) = a{x - xlnx) + {yi^i - yi.ilnyi^j) + (mi - mi In mi). 

i=l 

The coefficients (a, hi^i) are positive and defined from Ai as in the proof of Theorem 2.1 
from section 2.2. We also use the function Vee defined in (2.10) to consider 

n 

V{x, z) = To Vee{x, zi)+a^ VoFEizi) 

or equivalently 

n 

V{x, z) = To' Vee{x, z,)+aJ2 Pi ^* ( ^i.c | {-Ar') Zi ). 

Using the relation (2.19) and (2.11), we can compute the derivative of V along the 
trajectories of (3.2): 

n n 

^ = T'o'^ix, zi) + aTo ^ A Si ( Cicu | ) - o-Tq ^ Pix{eiu \ Zi) 

n 

+ a ^ /3i ( I ti^u ) (TqX- X*) . 

j=2 

Using Tq Xi = X* from the Definition 2.2 for the threshold we get 

n 

V = V$(x, Zi) + aJ2Pi{zi\ Bi,^ )x{T^- V) < 0. 

i=2 

By Liapunov theorem this ends the proof for the stability. The global asymptotic 
stability is obtained by a straightforward use of LaSalle's invariance principle, which 
ends the proof of Theorem 3.1. □ 

Remark 4. In the nongeneric case it can be shown, with the help of the Liapunov 
functions used in the theorem, that there exists a continuum of stable EE. We omit 
the proof. 

In the generic case, the dynamics of the system are completely determined. The 
nonnegative orthant is stratified in the union of stable manifolds corresponding to the 
different equilibria. Only the equilibrium corresponding to the winning strain has a 
basin of attraction with a nonempty interior. 

Remark 5. We have proved that the most virulent strain, that is, the strain 
which maximizes its respective threshold Tq, eliminates the other. We obtain the 
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same kind of result as in [7] , where the authors consider a SIR model with n strains 
of parasite. They consider that infection by one parasite strain excludes superinfection 
by other strains (this is also our case) and induces permanent immunity against all 
strains in case of recovery. They also guarantee limited population by considering a 
recruitment depending on the density in a monotone decreasing way. They find that 
the strain which maximizes the basic reproduction ratio eliminates the others. In the 
case considered by the authors, actually, using our notation, TZq = In fact in this 
model To and TZq coincide. This is also the case in our model when u = 0. Hence our 
result compares with the result of [7]. However in the case m ^ this is Tg , and not 
TZq, which distinguishes the fate of the strain. Our result is then different from [7], 
where this role is devoted to TZq. The same kind of remarks apply to [10] and [11]. 

Remark 6. In our model the chains are of equal length for each strain. If the 
chains are of unequal length, the proof is unchanged. We use equal length for no- 
tational convenience. A reason to have unequal length could be to model different 
behavior for two different strains of the parasite. 

4. Conclusion. In this article we have given a parasitic within-host model and 
have provided a stability analysis of this model. 

This model incorporates a number k of compartments for the parasitized target 
cells and considcirs n strains for the parasite. The rationale for including multicom- 
partments can be multiple. One reason is to take into account biological reasons, e.g., 
consideration of morphological or age classes. The second is for behavioral modeling 
reasons, e.g., to model delays described by gamma distribution functions. 

This model has been conceived from malaria infection, since it is well grounded 
that malaria is a multistrain infection. However other parasitic infections can be 
considered by this model. 

We prove that if the basic reproduction number satisfies T^o < 1, then the DFE 
is GAS; i.e., the parasite is cleared from the host. Our stability result when TZq > 1 
can be summarized as a competitive exclusion principle. To each i-strain we associate 
an individual threshold condition Tq as in Definition 2.2. If TZq > 1, if one strain has 
its individual threshold strictly larger than the thresholds of the other strains and if 
a mild sufficient condition is satisfied (for a constant recruitment, i.e., f{x) = A, this 
condition is simply u/3A < fJ-xl^'m), then there exists a GAS equilibrium on the 
positive orthant. This equilibrium corresponds to the extinction of all strains, except 
the strain with the largest threshold. This winning strain maximizes the threshold 
and not its individual basic reproduction number, which is different from previous 
analogous results of the literature. 
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